# Load necessary librarieslibrary(forecast)library(tseries)library(tidyverse)library(tidymodels)library(kableExtra)library(modeltime)library(timetk)library(cowplot)library(Metrics)# setting the seedset.seed(123)# Loading Data ----------------------------------------------------------------daily_deaths <-read_csv("data/new_daily_deaths.csv") %>%select(daily_deaths)east <-read_csv("data/east_daily.csv")midwest <-read_csv("data/midwest_daily.csv")south <-read_csv("data/south_daily.csv")west <-read_csv("data/west_daily.csv")
2. Data Splitting
Code
# Data Splitting ---------------------------------------------------------------east_splits <-time_series_split(east, initial ="1022 days", assess ="114 days")east_train <-training(east_splits)east_test <-testing(east_splits)midwest_splits <-time_series_split(midwest, initial ="1022 days", assess ="114 days")midwest_train <-training(midwest_splits)midwest_test <-testing(midwest_splits)south_splits <-time_series_split(south, initial ="1022 days", assess ="114 days")south_train <-training(south_splits)south_test <-testing(south_splits)west_splits <-time_series_split(west, initial ="1022 days", assess ="114 days")west_train <-training(west_splits)west_test <-testing(west_splits)## Convert Into Time Series ----------------------------------------------------east_train_ts <-ts(east_train %>%select(daily_deaths))east_test_ts <-ts(east_test %>%select(daily_deaths))east_ts <-ts(east %>%select(daily_deaths))midwest_train_ts <-ts(midwest_train %>%select(daily_deaths))midwest_test_ts <-ts(midwest_test %>%select(daily_deaths))midwest_ts <-ts(midwest %>%select(daily_deaths))south_train_ts <-ts(south_train %>%select(daily_deaths))south_test_ts <-ts(south_test %>%select(daily_deaths))south_ts <-ts(south %>%select(daily_deaths))west_train_ts <-ts(west_train %>%select(daily_deaths))west_test_ts <-ts(west_test %>%select(daily_deaths))west_ts <-ts(west %>%select(daily_deaths))## Plotting Training and Testing -----------------------------------------------par(mfrow=c(2,2))ggplot() +geom_line(data = east_test, aes(x = date, y = daily_deaths, color ="turquoise")) +geom_point(data = east_test, aes(x = date, y = daily_deaths, color ="turquoise"), size = .5) +geom_line(data = east_train, aes(x = date, y = daily_deaths, color ="orange")) +geom_point(data = east_train, aes(x = date, y = daily_deaths, color ="orange"), size = .5) +scale_color_manual(name="Data Set", values =c("turquoise", "orange"), labels =c("Training", "Testing")) +labs(y ="Daily Deaths",x ="Date",title ="East Region Training vs. Testing Data")
Code
ggplot() +geom_line(data = midwest_test, aes(x = date, y = daily_deaths, color ="turquoise")) +geom_point(data = midwest_test, aes(x = date, y = daily_deaths, color ="turquoise"), size = .5) +geom_line(data = midwest_train, aes(x = date, y = daily_deaths, color ="orange")) +geom_point(data = midwest_train, aes(x = date, y = daily_deaths, color ="orange"), size = .5) +scale_color_manual(name="Data Set", values =c("turquoise", "orange"), labels =c("Training", "Testing")) +labs(y ="Daily Deaths",x ="Date",title ="Midwest Region Training vs. Testing Data")
Code
ggplot() +geom_line(data = south_test, aes(x = date, y = daily_deaths, color ="turquoise")) +geom_point(data = south_test, aes(x = date, y = daily_deaths, color ="turquoise"), size = .5) +geom_line(data = south_train, aes(x = date, y = daily_deaths, color ="orange")) +geom_point(data = south_train, aes(x = date, y = daily_deaths, color ="orange"), size = .5) +scale_color_manual(name="Data Set", values =c("turquoise", "orange"), labels =c("Training", "Testing")) +labs(y ="Daily Deaths",x ="Date",title ="South Region Training vs. Testing Data")
Code
ggplot() +geom_line(data = west_test, aes(x = date, y = daily_deaths, color ="turquoise")) +geom_point(data = west_test, aes(x = date, y = daily_deaths, color ="turquoise"), size = .5) +geom_line(data = west_train, aes(x = date, y = daily_deaths, color ="orange")) +geom_point(data = west_train, aes(x = date, y = daily_deaths, color ="orange"), size = .5) +scale_color_manual(name="Data Set", values =c("turquoise", "orange"), labels =c("Training", "Testing")) +labs(y ="Daily Deaths",x ="Date",title ="West Region Training vs. Testing Data")
3. Testing for Stationarity
Code
### Using Augmented Dickey-Fuller Test -----------------------------------------adf_east <-adf.test(east_ts)print(adf_east) # p-value is 0.01 < 0.05, implying it is stationary
Augmented Dickey-Fuller Test
data: east_ts
Dickey-Fuller = -4.5566, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary
Code
adf_midwest <-adf.test(midwest_ts)print(adf_midwest) # p-value is 0.325 > 0.05, implying it is NOT stationary
Augmented Dickey-Fuller Test
data: midwest_ts
Dickey-Fuller = -2.5984, Lag order = 10, p-value = 0.325
alternative hypothesis: stationary
Code
diff_midwest_ts <-diff(midwest_ts)diff_adf_midwest <-adf.test(diff_midwest_ts)print(diff_adf_midwest) # p-value is 0.01 < 0.05, implying it is stationary with 1 diff
Augmented Dickey-Fuller Test
data: diff_midwest_ts
Dickey-Fuller = -13.43, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary
Code
adf_south <-adf.test(south_ts)print(adf_south) # p-value is 0.3018 > 0.05, implying it is NOT stationary
Augmented Dickey-Fuller Test
data: south_ts
Dickey-Fuller = -2.6533, Lag order = 10, p-value = 0.3018
alternative hypothesis: stationary
Code
diff_south_ts <-diff(south_ts)diff_adf_south <-adf.test(diff_south_ts)print(diff_adf_south) # p-value is 0.01 < 0.05, implying it is stationary with 1 diff
Augmented Dickey-Fuller Test
data: diff_south_ts
Dickey-Fuller = -11.36, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary
Code
adf_west <-adf.test(west_ts)print(adf_west) # p-value is 0.4098 > 0.05, implying it is NOT stationary
Augmented Dickey-Fuller Test
data: west_ts
Dickey-Fuller = -2.3981, Lag order = 10, p-value = 0.4098
alternative hypothesis: stationary
Code
diff_west_ts <-diff(west_ts)diff_adf_west <-adf.test(diff_west_ts)print(diff_adf_west) # p-value is 0.01 < 0.05, implying it is stationary with 1 diff
Augmented Dickey-Fuller Test
data: diff_west_ts
Dickey-Fuller = -12.247, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary
## ACF & PACF Plots ------------------------------------------------------------par(mfrow=c(3, 1))acf(east_ts, main ="East")pacf(east_ts, main ="East")plot(east_ts, main ="Daily Deaths: East")
Code
acf(midwest_ts, main ="Midwest")pacf(midwest_ts, main ="Midwest")plot(midwest_ts, main ="Daily Deaths: Midwest")
Code
acf(south_ts, main ="South")pacf(south_ts, main ="South")plot(south_ts, main ="Daily Deaths: South")
Code
acf(west_ts, main ="West")pacf(west_ts, main ="West")plot(south_ts, main ="Daily Deaths: South")
Code
plot(west_ts, main ="Daily Deaths: West")
5. Time Series Decomposition
Code
# total daily death decompositiondaily_ts <-ts(daily_deaths, frequency =365)decomposed_ts <-decompose(daily_ts)# Plot the original time series and the decomposed componentspar(mfrow=c(2,1))plot(daily_ts, main ="Total Original Time Series", ylab ="Value", col ="blue")plot(decomposed_ts$trend, main ="Trend Component", ylab ="Value", col ="red")
Code
plot(decomposed_ts$seasonal, main ="Seasonal Component", ylab ="Value", col ="green")plot(decomposed_ts$random, main ="Residual Component", ylab ="Value", col ="purple")
Code
# East Decomposition -----------------------------------------------------------east_daily <- east %>%select(daily_deaths)east_ts <-ts(east_daily, frequency =365)decomposed_east <-decompose(east_ts)# Plot the original time series and the decomposed componentspar(mfrow=c(2,1))plot(east_ts, main ="East Original Time Series", ylab ="Value", col ="blue")plot(decomposed_east$trend, main ="Trend Component", ylab ="Value", col ="red")
Code
plot(decomposed_east$seasonal, main ="Seasonal Component", ylab ="Value", col ="green")plot(decomposed_east$random, main ="Residual Component", ylab ="Value", col ="purple")
Code
## Midwest Decomposition ----------------------------------------------------------midwest_daily <- midwest %>%select(daily_deaths)midwest_ts <-ts(midwest_daily, frequency =365)decomposed_midwest <-decompose(midwest_ts)# Plot the original time series and the decomposed componentspar(mfrow=c(2,1))plot(midwest_ts, main ="Midwest Original Time Series", ylab ="Value", col ="blue")plot(decomposed_midwest$trend, main ="Trend Component", ylab ="Value", col ="red")
Code
plot(decomposed_midwest$seasonal, main ="Seasonal Component", ylab ="Value", col ="green")plot(decomposed_midwest$random, main ="Residual Component", ylab ="Value", col ="purple")
Code
## South Decomposition ----------------------------------------------------------south_daily <- south %>%select(daily_deaths)south_ts <-ts(south_daily, frequency =365)decomposed_south <-decompose(south_ts)# Plot the original time series and the decomposed componentspar(mfrow=c(2,1))plot(south_ts, main ="South Original Time Series", ylab ="Value", col ="blue")plot(decomposed_south$trend, main ="Trend Component", ylab ="Value", col ="red")
Code
plot(decomposed_south$seasonal, main ="Seasonal Component", ylab ="Value", col ="green")plot(decomposed_south$random, main ="Residual Component", ylab ="Value", col ="purple")
Code
## West Decomposition ----------------------------------------------------------west_daily <- west %>%select(daily_deaths)west_ts <-ts(west_daily, frequency =365)decomposed_west <-decompose(west_ts)# Plot the original time series and the decomposed componentspar(mfrow=c(2,1))plot(west_ts, main ="West Original Time Series", ylab ="Value", col ="blue")plot(decomposed_west$trend, main ="Trend Component", ylab ="Value", col ="red")
Code
plot(decomposed_west$seasonal, main ="Seasonal Component", ylab ="Value", col ="green")plot(decomposed_west$random, main ="Residual Component", ylab ="Value", col ="purple")
6. ARIMA Model
6.1 Initial Model
Code
# Building ARIMA Model ---------------------------------------------------------east_arima <-arima(east_train_ts, order=c(2,0,0))midwest_arima <-arima(midwest_train_ts, order=c(7,1,0))south_arima <-arima(south_train_ts, order=c(7,1,0))west_arima <-arima(west_train_ts, order=c(7,1,0))## Summary & Forecast ---------------------------------------------------------summary(east_arima)
Call:
arima(x = east_train_ts, order = c(2, 0, 0))
Coefficients:
ar1 ar2 intercept
0.8946 0.0181 198.4613
s.e. 0.0312 0.0313 39.1024
sigma^2 estimated as 12157: log likelihood = -6257.33, aic = 12522.66
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.02128508 110.2581 64.66103 -Inf Inf 0.9917366 -0.003359883
Code
summary(midwest_arima)
Call:
arima(x = midwest_train_ts, order = c(7, 1, 0))
Coefficients:
ar1 ar2 ar3 ar4 ar5 ar6 ar7
-0.6984 -0.6552 -0.6662 -0.6405 -0.625 -0.5942 0.2158
s.e. 0.0306 0.0326 0.0331 0.0336 0.033 0.0325 0.0305
sigma^2 estimated as 13720: log likelihood = -6315.71, aic = 12647.41
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.2590604 117.0771 67.37069 NaN Inf 0.4941415 -0.01395866
Code
summary(south_arima)
Call:
arima(x = south_train_ts, order = c(7, 1, 0))
Coefficients:
ar1 ar2 ar3 ar4 ar5 ar6 ar7
-0.7847 -0.7948 -0.782 -0.7725 -0.7874 -0.6062 0.0069
s.e. 0.0313 0.0349 0.035 0.0351 0.0349 0.0348 0.0312
sigma^2 estimated as 40704: log likelihood = -6869.67, aic = 13755.34
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.321387 201.654 117.0239 NaN Inf 0.5717169 -0.001802532
Code
summary(west_arima)
Call:
arima(x = west_train_ts, order = c(7, 1, 0))
Coefficients:
ar1 ar2 ar3 ar4 ar5 ar6 ar7
-0.7392 -0.7212 -0.7311 -0.7477 -0.7069 -0.5949 0.0674
s.e. 0.0312 0.0341 0.0343 0.0339 0.0342 0.0340 0.0311
sigma^2 estimated as 8677: log likelihood = -6080.76, aic = 12177.52
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.2961085 93.10308 54.88534 -Inf Inf 0.5462533 -0.006060256
Code
east_values <-forecast(east_arima, h =length(east_test_ts))midwest_values <-forecast(midwest_arima, h =length(midwest_test_ts))south_values <-forecast(south_arima, h =length(south_test_ts))west_values <-forecast(west_arima, h =length(west_test_ts))
6.2 Accuracy
Code
## Set Up Plotseast_test_plot <-as.data.frame(east_test_ts) %>%mutate(num =seq(1023, (1023+113)))midwest_test_plot <-as.data.frame(midwest_test_ts) %>%mutate(num =seq(1023, (1023+113)))south_test_plot <-as.data.frame(south_test_ts) %>%mutate(num =seq(1023, (1023+113)))west_test_plot <-as.data.frame(west_test_ts) %>%mutate(num =seq(1023, (1023+113)))
### South Regionsouth_initial_arima_acc <- forecast::accuracy(south_values)south_initial_arima_acc %>%kbl(caption ="<center>Initial South Accuracy<center>") %>%kable_styling()
### West Regionwest_initial_arima_acc <- forecast::accuracy(west_values)# west_initial_arima_acc %>% # kbl(caption = "<center>Initial West Accuracy<center>") %>% # kable_styling()west_pred <-as.data.frame(west_values$mean)west_pred2 <-cbind(west_pred, west_test_plot)west_initial_arima_mae <-mae(west_pred2$daily_deaths, west_pred2$x)west_initial_arima_mase <-mase(west_pred2$daily_deaths, west_pred2$x)west_initial_arima_plot <-ggplot(west_pred2) +geom_line(aes(x = num, y = daily_deaths), color ="dodgerblue", linewidth = .75) +geom_line(aes(x = num, y = x), color ="orange", linewidth = .75) +scale_x_continuous(limits =c(1023,1136),breaks =c(1023, 1069, 1115),labels =c("Nov 2022", "Jan 2023", "March 2023")) +labs(title ="West Region Initial ARIMA",y =NULL,x =NULL) +scale_color_manual(values =c("orange", "dodgerblue"), labels =c("Actual", "Predicted"))+theme_bw()west_initial_arima_plot
6.3 Grid Search
Code
# Grid Search for p and q ------------------------------------------------------p_values <-0:8# AR orderd_values <-0:0# I orderq_values <-0:8# MA order## East - Perform grid search for ARIMA parameters ----------------------------best_east_model <-NULLbest_east_aic <-Inffor (p in p_values) {for (d in d_values) {for (q in q_values) {# Fit ARIMA model arima_model <-tryCatch( { fit <-Arima(east_ts, order =c(p, d, q)) fit },error =function(e) {NULL } )# Check if ARIMA model was successfully fittedif (!is.null(arima_model)) {# Check AIC criterionif (AIC(arima_model) < best_east_aic) { best_east_model <- arima_model best_east_aic <-AIC(arima_model) best_east_order <-c(p, d, q) } } } }}
Code
## Midwest - Perform grid search for ARIMA parameters -------------------------p_values <-0:8# AR orderd_values <-1# I orderq_values <-0:8# MA orderbest_midwest_model <-NULLbest_midwest_aic <-Inffor (p in p_values) {for (d in d_values) {for (q in q_values) {# Fit ARIMA model arima_model <-tryCatch( { fit <-Arima(midwest_ts, order =c(p, d, q)) fit },error =function(e) {NULL } )# Check if ARIMA model was successfully fittedif (!is.null(arima_model)) {# Check AIC criterionif (AIC(arima_model) < best_midwest_aic) { best_midwest_model <- arima_model best_midwest_aic <-AIC(arima_model) best_midwest_order <-c(p, d, q) } } } }}
Code
## South - Perform grid search for ARIMA parameters ---------------------------p_values <-0:8# AR orderd_values <-1# I orderq_values <-0:8# MA orderbest_south_model <-NULLbest_south_aic <-Inffor (p in p_values) {for (d in d_values) {for (q in q_values) {# Fit ARIMA model arima_model <-tryCatch( { fit <-Arima(south_ts, order =c(p, d, q)) fit },error =function(e) {NULL } )# Check if ARIMA model was successfully fittedif (!is.null(arima_model)) {# Check AIC criterionif (AIC(arima_model) < best_south_aic) { best_south_model <- arima_model best_south_aic <-AIC(arima_model) best_south_order <-c(p, d, q) } } } }}
Code
## West - Perform grid search for ARIMA parameters ----------------------------p_values <-0:8# AR orderd_values <-1# I orderq_values <-0:8# MA orderbest_west_model <-NULLbest_west_aic <-Inffor (p in p_values) {for (d in d_values) {for (q in q_values) {# Fit ARIMA model arima_model <-tryCatch( { fit <-Arima(west_ts, order =c(p, d, q)) fit },error =function(e) {NULL } )# Check if ARIMA model was successfully fittedif (!is.null(arima_model)) {# Check AIC criterionif (AIC(arima_model) < best_west_aic) { best_west_model <- arima_model best_west_aic <-AIC(arima_model) best_west_order <-c(p, d, q) } } } }}
# Building SARIMA model ------------------------------------------------------east_sarima <-arima(east_train_ts, order =c(7,0,7), seasonal =list(order =c(1,0,1)))midwest_sarima <-arima(midwest_train_ts, order =c(8,1,8), seasonal =list(order =c(1,1,1)))south_sarima <-arima(south_train_ts, order =c(6,1,7), seasonal =list(order =c(1,1,1)))west_sarima <-arima(west_train_ts, order =c(8,1,8), seasonal =list(order =c(1,1,1)))# Summary of the SARIMA modelsummary(east_sarima)
## Forecasting with the SARIMA model -------------------------------------------east_values <-forecast(east_sarima, h =length(east_test_ts))midwest_values <-forecast(midwest_sarima, h =length(midwest_test_ts))south_values <-forecast(south_sarima, h =length(south_test_ts))west_values <-forecast(west_sarima, h =length(west_test_ts))
7.2 Accuracy
Code
## Plotting & Accuracy## East Regioneast_initial_sarima_acc <- forecast::accuracy(east_values)# east_initial_sarima_acc %>% # kbl(caption = "<center>Initial SARIMA East Accuracy<center>") %>% # kable_styling()east_pred <-as.data.frame(east_values$mean)east_pred2 <-cbind(east_pred, east_test_plot)east_initial_sarima_mae <-mae(east_pred2$daily_deaths, east_pred2$x)east_initial_sarima_mase <-mase(east_pred2$daily_deaths, east_pred2$x)### East Regioneast_initial_sarima_plot <-ggplot(east_pred2) +geom_line(aes(x = num, y = daily_deaths, color ="dodgerblue"), linewidth = .75) +geom_line(aes(x = num, y = x, color ="orange"), linewidth = .75) +scale_x_continuous(limits =c(1023,1136),breaks =c(1023, 1069, 1115),labels =c("Nov 2022", "Jan 2023", "March 2023")) +labs(title ="East Region Initial SARIMA",y =NULL,x =NULL) +scale_color_manual(values =c("dodgerblue", "orange"), labels =c("Actual", "Predicted"))+theme_bw() +theme(legend.title =element_blank())east_initial_sarima_plot
## Define the parameter grids --------------------------------------------------p_grid <-0:6# AR parameterd_grid <-0:0# Differencing parameterq_grid <-0:6# MA parameterP_grid <-0:1# Seasonal AR parameterD_grid <-0:0# Seasonal differencing parameterQ_grid <-0:1# Seasonal MA parameters <-7# Seasonal period## East Grid Search ------------------------------------------------------------# Create East Regioneast_train <- east_train %>%select(daily_deaths)# Create an empty data frame to store the resultseast_results <-data.frame(order =character(),seasonal =character(),AIC =numeric(),stringsAsFactors =FALSE)# Perform grid searchfor (p in p_grid) {for (d in d_grid) {for (q in q_grid) {for (P in P_grid) {for (D in D_grid) {for (Q in Q_grid) { order <-c(p, d, q) seasonal <-c(P, D, Q, s) model <-tryCatch( { fit <-arima(east_train, order = order, seasonal = seasonal) AIC_val <-AIC(fit)# Append the results to the data frame east_results <-rbind(east_results, data.frame(order =paste(order, collapse =", "),seasonal =paste(seasonal, collapse =", "),AIC = AIC_val)) },error =function(e) NULL ) } } } } }}d_grid <-1# Differencing parameterD_grid <-1# Differencing parameter## Midwest Grid Search ------------------------------------------------------------# Create Midwest Regionmidwest_train <- midwest_train %>%select(daily_deaths)# Create an empty data frame to store the resultsmidwest_results <-data.frame(order =character(),seasonal =character(),AIC =numeric(),stringsAsFactors =FALSE)# Perform grid searchfor (p in p_grid) {for (d in d_grid) {for (q in q_grid) {for (P in P_grid) {for (D in D_grid) {for (Q in Q_grid) { order <-c(p, d, q) seasonal <-c(P, D, Q, s) model <-tryCatch( { fit <-arima(midwest_train, order = order, seasonal = seasonal) AIC_val <-AIC(fit)# Append the results to the data frame midwest_results <-rbind(midwest_results, data.frame(order =paste(order, collapse =", "),seasonal =paste(seasonal, collapse =", "),AIC = AIC_val)) },error =function(e) NULL ) } } } } }}## South Grid Search ------------------------------------------------------------# Create South Regionsouth_train <- south_train %>%select(daily_deaths)# Create an empty data frame to store the resultssouth_results <-data.frame(order =character(),seasonal =character(),AIC =numeric(),stringsAsFactors =FALSE)# Perform grid searchfor (p in p_grid) {for (d in d_grid) {for (q in q_grid) {for (P in P_grid) {for (D in D_grid) {for (Q in Q_grid) { order <-c(p, d, q) seasonal <-c(P, D, Q, s) model <-tryCatch( { fit <-arima(south_train, order = order, seasonal = seasonal) AIC_val <-AIC(fit)# Append the results to the data frame south_results <-rbind(south_results, data.frame(order =paste(order, collapse =", "),seasonal =paste(seasonal, collapse =", "),AIC = AIC_val)) },error =function(e) NULL ) } } } } }}## West Grid Search ------------------------------------------------------------# Create West Regionwest_train <- west_train %>%select(daily_deaths)# Create an empty data frame to store the resultswest_results <-data.frame(order =character(),seasonal =character(),AIC =numeric(),stringsAsFactors =FALSE)# Perform grid searchfor (p in p_grid) {for (d in d_grid) {for (q in q_grid) {for (P in P_grid) {for (D in D_grid) {for (Q in Q_grid) { order <-c(p, d, q) seasonal <-c(P, D, Q, s) model <-tryCatch( { fit <-arima(west_train, order = order, seasonal = seasonal) AIC_val <-AIC(fit)# Append the results to the data frame west_results <-rbind(west_results, data.frame(order =paste(order, collapse =", "),seasonal =paste(seasonal, collapse =", "),AIC = AIC_val)) },error =function(e) NULL ) } } } } }}
# Rebuilding SARIMA Models -----------------------------------------------------east_sarima <-arima(east_train_ts, order =c(4, 0, 6), seasonal =list(order =c(1, 0, 0), period =7))midwest_sarima <-arima(midwest_train_ts, order =c(3, 1, 6),seasonal =list(order =c(0, 1, 1), period =7))south_sarima <-arima(south_train_ts, order =c(5, 1, 6),seasonal =list(order =c(1, 1, 0), period =7))west_sarima <-arima(west_train_ts, order =c(6, 1, 5),seasonal =list(order =c(0, 1, 1), period =7))## Forecasting with the NEW SARIMA model ---------------------------------------east_values <-forecast(east_sarima, h =length(east_test_ts))midwest_values <-forecast(midwest_sarima, h =length(midwest_test_ts))south_values <-forecast(south_sarima, h =length(south_test_ts))west_values <-forecast(west_sarima, h =length(west_test_ts))
Code
## Plotting & Accuracy## East Regioneast_tuned_sarima_acc <- forecast::accuracy(east_values)# east_tuned_sarima_acc %>% # kbl(caption = "<center>Tuned SARIMA East Accuracy<center>") %>% # kable_styling()east_pred <-as.data.frame(east_values$mean)east_pred2 <-cbind(east_pred, east_test_plot)east_tuned_sarima_mae <-mae(east_pred2$daily_deaths, east_pred2$x)east_tuned_sarima_mase <-mase(east_pred2$daily_deaths, east_pred2$x)### East Regioneast_tuned_sarima_plot <-ggplot(east_pred2) +geom_line(aes(x = num, y = daily_deaths, color ="dodgerblue"), linewidth = .75) +geom_line(aes(x = num, y = x, color ="orange"), linewidth = .75) +scale_x_continuous(limits =c(1023,1136),breaks =c(1023, 1069, 1115),labels =c("Nov 2022", "Jan 2023", "March 2023")) +labs(title ="East Region Tuned SARIMA",y =NULL,x =NULL) +scale_color_manual(values =c("dodgerblue", "orange"), labels =c("Actual", "Predicted"))+theme_bw() +theme(legend.title =element_blank())east_tuned_sarima_plot
# Building AutoARIMA Model -----------------------------------------------------east_arima <-auto.arima(east_train_ts, seasonal =TRUE)midwest_arima <-auto.arima(midwest_train_ts, seasonal =TRUE)south_arima <-auto.arima(south_train_ts, seasonal =TRUE)west_arima <-auto.arima(west_train_ts, seasonal =TRUE)# Summary of the ARIMA modelsummary(east_arima)
Series: east_train_ts
ARIMA(5,1,2)
Coefficients:
ar1 ar2 ar3 ar4 ar5 ma1 ma2
0.5142 -0.7304 -0.1394 -0.3479 -0.2207 -1.0295 0.7353
s.e. 0.0480 0.0344 0.0426 0.0329 0.0418 0.0410 0.0265
sigma^2 = 6861: log likelihood = -5956.39
AIC=11928.79 AICc=11928.93 BIC=11968.22
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1506172 82.50546 47.1043 NaN Inf 0.7224607 0.01707355
Code
summary(midwest_arima)
Series: midwest_train_ts
ARIMA(5,1,3)
Coefficients:
ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3
0.1952 -0.6269 -0.3389 -0.2999 -0.4563 -1.2016 0.9332 -0.3022
s.e. 0.0607 0.0312 0.0356 0.0276 0.0442 0.0702 0.0718 0.0382
sigma^2 = 21712: log likelihood = -6544.65
AIC=13107.31 AICc=13107.49 BIC=13151.67
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.3932594 146.6988 85.10252 NaN Inf 0.6241985 -0.03369519
Code
summary(south_arima)
Series: south_train_ts
ARIMA(5,1,4)
Coefficients:
ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3
0.1881 -0.9049 -0.0800 -0.4542 -0.5567 -0.9401 1.0118 -0.7340
s.e. 0.0349 0.0292 0.0464 0.0251 0.0303 0.0357 0.0434 0.0525
ma4
0.4897
s.e. 0.0402
sigma^2 = 40663: log likelihood = -6864.82
AIC=13749.64 AICc=13749.86 BIC=13798.93
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1420877 200.6629 119.7558 NaN Inf 0.5850634 -0.06523796
Code
summary(west_arima)
Series: west_train_ts
ARIMA(4,1,3)
Coefficients:
ar1 ar2 ar3 ar4 ma1 ma2 ma3
-0.0876 0.1258 -0.5602 -0.4779 -0.5985 -0.3463 0.5343
s.e. 0.0505 0.0450 0.0330 0.0309 0.0509 0.0784 0.0421
sigma^2 = 11991: log likelihood = -6241.45
AIC=12498.9 AICc=12499.04 BIC=12538.33
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1543984 109.0727 72.43011 NaN Inf 0.7208697 -0.08454213
Code
## Forecasting Test Data -------------------------------------------------------east_values <-forecast(east_arima, h =length(east_test_ts))midwest_values <-forecast(midwest_arima, h =length(midwest_test_ts))south_values <-forecast(south_arima, h =length(south_test_ts))west_values <-forecast(west_arima, h =length(west_test_ts))